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Abstract. Simulations of standard ID and 2D quantum walks have 
been performed within Quantum Computer Simulator (QCS system) en- 
vironment and with the use of GPU supported by CUDA technology. In 
particular, simulations of quantum walks may be seen as an appropriate 
benchmarks for testing calculational power of the processors used. It was 
demonstrated by a series of tests that the use of CUDA based technology 
radically increases the computational power compared to the standard 
CPU based computations. 
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1 Introduction 

Recently the concept of quantum walks attracted a big attention as they pro- 
vide us with a very promising source of ideas for constructing new quantum 
algorithms |15j . In particular exponential speedups of some classical problems 
have been discovered again (like in the Shor's algorithm case), the exponen- 
tially faster hitting described in [3] and [5] seems to be a good example for 
this. Additionally, certain although not so particular as exponential speedups 
are, speedups of some classical problems like k-distinction problem, triangle and 
k-cliquc algorithms are examples of them [5] , [T5] , [Hj have been obtained by the 
use of quantum walks concept. Another inspirations for studying quantum walks 
based applications may come from the observation [5] that quantum walks might 
play the role of universal quantum calculations tool. 

It is the main purpose of the present contribution to test some particular 
properties of quantum walks by using certain simulating tools. The main tools 
used for our simulations are the Zielona Gora Quantum Computer Simulator (as 
the main tool) and the GPGPU equipped with new computational technology 
offered by CUDA (Compute Unified Device Architecture). 

It is one of the main result of this note the demonstration of powerful incre- 
mentation of calculational abilities if we use GPGPU of new generation compared 
with standard CPU based computations. 

Organisation of this note is as follows: in Section ([2]) we provide a reader 
with the basic definitions and constructions connected with quantum walks on 



general graphs. Numerical implementations of quantum walks on certain lattice 
structures coming from that of Z 1 and T? together with the corresponding nu- 
merical algorithms are presented in Section Additionally several examples 
of our simulations will be presented and discussed briefly there. 

2 Some mathematical preliminaries 

Let G = (V, E, I) be a graph and where as obviously V stands for the set of 
vertices of G the number of elements |V| of which is equal to N, E stands for 
the set of (undirected) edges of G and I : E — > {V, V} is the edges labelling 
function. The corresponding incidence matrix of G will denoted as Mq. With 
the use of Mq the continuous time Markov walk on G can be described by 
the corresponding, canonical Markovian transition semi-group T t = exp(— tMa) 
and (eventually) an initial distribution ttq. For any v G V we denote by d(v) 
the corresponding degree of vertex v and let then P(«) = (Pi,P 2 , . . . , Pd(v)) be 
discrete probability measure assigning a probability that the walker jump with 
probability pt by the use of i-th edge e$ connecting the vertex v with 1(a)(2) = u. 
The system (P, no), where ttq is an initial distribution gives rise to the discrete 
step random walk on G. 

The corresponding quantum walks on G can be constructed in the following 
way. 

The continuous time quantum walk on G: by the very definition starting from 
the vertex vq (with probability ttq) after time t we arrive at the vertex obtained 
by the action of the unitary group Uf = exp(itHc) (where Hq stands for the 
corresponding graph Hamiltonian) acting in the Hilbert space VP = ® v£ y\v). 
However only the discrete time processes will be discussed in this note. 

2.1 The discrete time Markovian quantum walks 

LctH G = ®vev\v) be the canonical Hilbert space associated with G, obviously 
H G ~ C N . For any v <E V let d(v) be degree of v. Then the local Hilbert space 
H v is defined as a space isomorphic to C d ( v \ explicitly H v = © e Je u ), e u runs 
over all edges connecting the vertex v with others. A collections C = (C v ,v G 
V) of unitary maps acting on the spaces % v , v G V and fulfilling additionally 
certain natural coincidence conditions, see i.e. [B]; will be called a "coin flip 
transformation sequence" . In other words, for any v G V: 



The global Hilbert space H is defined as TL = ® v ev\v) ® "Hv and the corre- 
sponding discrete quantum walk on G, providing the family C is given, can be 
defined by its one step transformation: 



C v : \v)®H v -> |u> ®Uv 
where \v) ® \u) — > \v) <g> C v \uj) 



(1) 



U = S(I®C) 



(2) 



where the shift transformation S is defined as: 

S\v,e) = \v',e) if 1(e) = {«,«'}. (3) 

Several questions (with analogy to the classical case, especially the prob- 
lems connected to the mixing and hitting times on a large class of graphs have 
been studied intensively) can be studied, the question about limiting probability 
distributions and hitting times are the most popular one [5]. 

Although intensive simulations of quantum walks on many complex graphs 
are planned to be done we have concentrated first on some simplest quantum 
walks on the infinite (the finite amount of with appropriate boundary condition 
are considered in real simulation tasks of course) graphs Z 1 and Z 2 that we 
describe now. 



2.2 Quantum walks on lattice Z 1 and on lattice Z 2 

With the lattice Z 1 we associate the Hilbert space /2(Z) = ® n £i\n) and the coin 
flip transformation C acting in C 2 = \R) © \L) symbolising the possible steps in 
the right (\R)) or left direction (\L)) is given. Then the corresponding Hilbert 
space /2(Z) ® C 2 can be seen as a space of infinite sequences \tp) « ((ay, k), oij 6 
Z, k G {R, L}), i.e. any vector \t/j) £ ^(Z) ® C 2 can be given by: 

= a jk\j,k) where ^|ay ife | 2 = l (4) 

jeZ,k=L,R j.k 

Different choices of C and shift operators lead to different models of quantum 
walks on the lattice Z 1 . 

The graph Hilbert space Tl z for the 2D lattice Z 2 is defined as 

R U2 = (B J:m ez\j,m) = Z 2 (Z 2 ). 

The degrees of all vertices are equal to 4 and therefore the local Hilbert spaces 
are isomorphic with C 4 to be identified with R,L,U,D (right, left, up, down) 
steps on the lattice. The total space H = ^(Z 2 ) ® C 4 and the typical vector 
|?A) £ % can be decomposed as 

4 +oo 4 +oo 

1^) = ^2 ^2 Uj,k,rn,n\j,k)\m,n) with ^ ^2 \ a 3,k,m,n\ 2 = 1 (5) 
j : k—l m.n— — oo j.k—1 m.n— — oo 

A different versions of the corresponding coin flip transformation C and the 
shift operators (reflecting some additional topological constraints) then lead to 
different quantum walk models on Z 2 lattice. Some of them will be presented for 
simulations performed in the next section including some quantum walk models 
on Z 1 as well. 



3 The algorithm for simulating quantum walks on 
GPGPU 



The calculation routine for simulation of quantum walks can be build directly 
from the definition of state of quantum walker walking on the lattice Z 2 . In 
general a state of the quantum walker at time t is given in the following way: 

1 +00 

IVK*)) = X X a j,k,m, n (t)\j,k)\m,n) ( 6 ) 

j,k=0 m,n— — 00 

where a^ k , m ,n{t) 6 C and J2j.kJ2 m ,n \ a j,k,m,n{t)\ = 1- The evolution of the 
quantum walker system over time t is expressed by following unitary operator 

U = S{C®I) (7) 

where S is the shift operator, / represents the identity operator and C is the coin 
operator (in most cases we can assume that the coin is represented by Hadamard 
operator, but there exist other representations of the coin operator e.g. Fourier 
and Grover coins) which acts on the local H2 65 H2 subspace of whole walker 
Hilbcrt space system ^(Z 2 ) ® H2 ® H2 ■ 

In this contribution we propose rather special definition of shift operator for 
two-dimensional quantum walks. A comparison of our definition with those used 
frequently can be found in |10j and [8]. We also use the random broken links 
(termed RBL) technique, first developed for one dimensional quantum walks and 
introduced in [TJ] . The RBL technique was generalised for two-dimensional case 
in [TO]. 

The used definition of shift operator which coincides with physical and math- 
ematical lattice is the following 

1 +00 

5= E X \3,k}(j,k\®\m + (-iy(l-6 jik ),n + (-iy5 jik )(m,n\ (8) 

j,k— m,n— — 00 

where is Dirac discrete delta function. 

The following function includes the possibility of appearance of the broken 
line in a path in between site (j, k) and (m, n) in example depicted on Fig. Q] 

f (-l) j - if link to m + - Sj, k ), 

TZBC(j, k, m,n) = < n + (— l) 3 '<5j,k is closed (9) 

[0 — if link is open 

where j, k e {0, 1}. 

After applying shift operator ((5J to state © the evolution can be summarised 
in following way: 

1 

+TZBC(j,k,m,n),k(BTZBC(j,k,m,,n),j' ,k' 

j',k'=0 

(*) (io) 
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Fig. 1. The idea of broken links for one-dimensional quantum walk (a) where 
we show a broken links between sites, part (b) and part (c) represent examples 
of natural and diagonal lattices for two-dimensional quantum walk 



where © represents addition modulo two. 

The algorithm to simulate two-dimensional quantum walks is directly basing 
on the evolution given by equation (|10[) . It can be implemented on the traditional 
architecture for standard CPU and of course for GPGPU based solutions. 

Figure © shows a fragment of calculation routine for GPGPU which task is 
to compute values for the next iteration of quantum walk in the segment case. 
In each iteration all points in the segment attain a new values. This means that 
the all available GPGPU cores are fully used. Additionally, the efficiency can be 
increased by better usage of threads available in CUDA architecture. 



_global void one_iteration_segment ( 

cuFloatComplex *A0, cuFloatComplex *A1, 

cuFloatComplex *AtempO, cuFloatComplex *Atempl, int *RBL0, 

int *RBL1, cuFloatComplex *C, int N) { 
int m = blockldx.x * blockDim.x + threadldx.x; 
int L, cidxl, cidx2; 

if (m<N) { 

L = RBLO [m] ; cidxl=(L*2); cidx2=(L*2)+l ; 

Atempl[m] = cuCaddf (cuCaddf (Atempl [m] , cuCmulf (C [cidxl] , A0[m+L])), 

cuCaddf (Atempl [m] , cuCmulf (C [cidx2] , Al [m+L] ) ) ) : 

L = RBL1 [m] ; cidxl=( (1+L) *2) ; cidx2=( (1+L) *2) + l ; 

AtempO[m] = cuCaddf (cuCaddf (At empO [m] , cuCmulf (C [cidxl] , A0[m+L])), 

cuCaddf (AtempO [m] , cuCmulf (C [cidx2] , Al [m+L] ) ) ) : 



Fig. 2. The kernel function to compute trajectories of one-dimensional quantum 
walks on segment 



The function to calculate the probability distribution for one-dimensional 
quantum walk on the line is very similar to the segment case with one important 
difference. In the i-th iteration the quantum walker cannot be farther than i sites 
from its initial position. The necessary change in GPGPU routine is expressed 
as 

int m = blockldx.x * blockDim.x + threadldx.x; 
int left , right ; 

left = max (midpoint - extra - iteration, 1); 
right = min (midpoint + extra + iteration, N-l) ; 
if (m>=left && m<=right) { . . . } 

The probability distributions connected to a calculated trajectory for segment 
and line are not the same, which is illustrated on the Fig. 
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Fig. 3. The one-dimensional quantum walk on segment (left plot) and on line 
(right plot). In both cases amount of iterations are the same and equal to 100 




Fig. 4. Probability distribution of two-dimensional quantum walks with Fourier 
walker and Fourier's initial state on the diagonal lattice depicted on plot (a) and 
natural lattice depicted on plot (b) 



3.1 Complexity analysis briefly 



The computational complexity T ld of simulation of quantum walk for one- 
dimensional case strictly depends on the size of segment or line used. Let TV 
be a size of segment for arising in i-th iteration, then: 

iV-l i / 1 \ 

m=0j=0 V fe=0 / 

It is possible to write equation (jllj) in such form because operations RBL m and 
OP m need the same constant time of work which is independent on the number 
of iteration. The symbol T&b L m stands for amount of time necessary to process 
information about broken links, the second symbol Top m stands for amount of 
time necessary to process the probability amplitude of selected point. The use 
of GPGPU means, that the processing time can be divided by the number N c 
of available cores, because operations RBL m and OP m are independent for each 
point in the segment or line (as well as in two-dimensional case): 

T^pgpuM = Tl d (N)/N c (12) 

The simulation of quantum walk on line shows one important difference com- 
paring to segment case. The point is that the new values of points are calculated 
only in the partially and the size of simulated portion strictly depends on the 
number of iteration. The complexity can be denoted as (where TpBL m and Top m 
have the same sense as in previous segment case): 

T? d (N) =J2I1 ( T RBL m Y ( T OP m ) \ =4-T RBLm -To Pm -((n-l i ) + l) (13) 

m=lj j=0 V fc=0 / 

where the U and are given by: 

k = max(MP - ES - i, 1), r< = min(MP + ES + i, N - 1) (14) 

where MP is the midpoint index of line, the value ES is used to pad and therefore 
prevents from range errors. However, these values are constant so the difference 
of Ti and U for i-th iteration can be expressed as 

(n - k) = (2 • i) + 1. (15) 

In the case of small systems (ten, twenty or fifty points), the equation (fl~5j) shows 
that the most of available cores in GPGPU are not fully utilised. In the case of 
bigger systems this problem is not arising. The difference in (|15|) for higher 
iteration number is bigger and what is more the values of this difference for all 
iterations form the arithmetic sequence, what means that in sense of complexity 
theory that only linear speedup is achieved, however for segment with size of 
10000 points the obtained speedup is nearly hundredfold. The corresponding 
times have been depicted in Tables ([1]), ^ and 



Table 1. The measured times of calculations of one-dimensional quantum walks 
trajectories for segments with different sizes (without broken links) 





Core 2 Duo 8400 (1 core) 


Geforce 9600 GT (64 cores) 


Geforce 280 (240 cores) 


Size 


Time in ms 


Time in ms 


Time in ms 


100 


6.00 


4.51 


2.48 


1000 


663.00 


40.186 


24.171 


5000 


22685.00 


612.67 


281.57 


10000 


96362.00 


2324.3002 


884.297 



Table 2. The measured times of calculations of two-dimensional quantum walks 
trajectories for diagonal lattice without broken links. The measured times for 
the case of Core 2 Duo and two-thread computational routine are presented in 
bracket 





Core 2 Duo 8400 (1 core) 


Geforce 9600 GT (64 cores) 


Geforce 280 (240 cores) 


Size 


Time in ms 


Time in ms 


Time in ms 


100 


700 (413) 


195 


64 


200 


5480 (2945) 


1003 


363 


300 


19140 (10512) 


3423 


1137 


400 


47230 (25785) 


8123 


3127 


500 


92530 (50274) 


14462 


4706 



Table 3. The measured times of the calculations of two-dimensional quantum 
walks trajectories for diagonal lattice without broken links. The simulation was 
performed on two Intel Xeon E5420 2.50 Ghz processors, the multi-threaded 
calculation subroutine was compiled with GCC compiler with "-03" option 





(1-thread) 


(2th) 


(4th) 


(8th) 


Size 


Time in ms 


Time in ms 


Time in ms 


Time in ms 


100 


770 


409 


299 


206 


200 


6740 


3281 


2190 


1382 


300 


22910 


12029 


7130 


4848 


400 


56610 


29853 


17699 


12069 


500 


114060 


65292 


33275 


24602 



The use of equation (| 10[) allows to estimate computational complexity of 
two-dimensional quantum walk trajectory calculations in i-th iteration: 

m=\hi n=\bi j=0 k=0 

16 • T RBL(m n j k) ■ T OP(mMk) • ((rbi - lb«) + l) 2 , (16) 

where N means the length of trajectory calculated. 

The variables lb^ and rb^ have the same meaning as k, Ti introduced before 
and are calculated in the following way: 

lbi = max(MP-ES-i,l), rbi = min(MP + ES + i, 2 ■ MP - 1). (17) 

The expression ((rbi — lbj) + l) 2 in equation (|16|) can be expressed in the 
following way 

((rb, - lb 4 ) + l) 2 = ((2-i) + l) 2 . (18) 

Spccdups obtained for 2D quantum walks are presented in Fig. arc based on 
results presented in Tables (HJ and 
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Fig. 5. The obtained values of speedup for simulations of two dimensional quan- 
tum walks on diagonal lattice. Graphs labelled as 2th, 4th and 8th are compared 
to one-thread computational routine. Graphs marked as 64th and 240th repre- 
sent the speedup obtained by computational routine executed on the Geforce 
video card 9600 GT and GTX 280 respectively 



4 Conclusions and further work 

The specialised software to simulate one and two dimensional random quantum 
walks without and with broken links has been presented in this article. The used 
software is a part of the Quantum Computing Simulator presented in |13j . A 
significant speedup of the simulations process comparing to previous paper [8] 



have been achieved. The used technologies enables to simulate effectively much 
more complex quantum walks then previously known. Additionally, certain more 
deeper notions connected to the analysis of quantum walks behaviour can be 
analysed using computer simulations as an appropriate tool. 
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